% Subroutine to select band 4 (GMS05) time series
% subsection for processing
% Uses artificial neural networks with saved weigths
% First NN checks the pattern quality and
% the second one gives overall quality
% Ref: C. Manoj & N. Nagarajan,2003. Geophys. J. Intl. 153, 1-15.
% Latest date 06.03.2003

function[ProcDef] = SelANN(TS,ProcDef);

fprintf('SelANN -> ANN editing\n');
load c:\manoj\data\ann\bpns001.mat;    %--weight file(s)
load c:\manoj\data\ann\bpnStmp.mat; 
matrix = TS.matrix;
clear TS;
trf=['logsig';'logsig'];
%---------- declaration -------------------------------
block = 256;
TAP10 = [0.9877 0.9511 0.8910 0.8090 0.7071 0.5878 0.4540 0.3090 0.1564 0.0000];
TAP20 = [0.0000 0.1564 0.3090 0.4540 0.5878 0.7071 0.8090 0.8910 0.9511 0.9877];
ntaper = 10;
Nsq = block*block;
AA = size(matrix);
length = AA(1)/block;
start = 0;
st = start*block + 1;
en = (start+1)*block;
nchannel = 5;
choice = 'ry'; 
Corr1= [1 4]; % Channels for correlation
Corr2= [2 3]; % ex = 1 ey = 2 hx = 3 hy = 4 
AB = size(W1);
len = AB(2);    % number of harmonics for bpn fwd cal
selection = [0 0.9 0.1 0];
fact = 0.0;

%---------- actual stuff ------------------------------



for i = (start+1):length,

for u = 1:5,
        K = matrix(st:en,u);
	bias  = mean(K);
	P2 = sum(K(1:block/2));
	P1 = sum(K(1+block/2:block));
	trend = 4*(P1-P2)/Nsq;

%---------Trend & bias removal-------------------------
	for n=1:block,
		K(n) = K(n) - (bias +trend*n);
	end;
	MI = min(K);
	MA = max(K);
	AMP(u) = MA-MI;

%---------Normalization--------------------------------
	for n=1:block,
		K(n) = (K(n) - MI)/(MA - MI);
		end;
	Channel(u,:) = K';
	end; 	
%---------------tapering the edges--------------------------
MS =  (Channel(1:nchannel,1:ntaper)');
ME =  (Channel(1:nchannel,block-ntaper+1:block)');
for ii = 1:nchannel,
tapmat(ii,:) = (MS(:,ii).*TAP20')';
end;
Channel(1:nchannel,1:ntaper) 		=	tapmat;
for ii = 1:nchannel,
tapmat(ii,:) = (ME(:,ii).*TAP10')';
end;
Channel(1:nchannel,block-ntaper+1:block) = 	tapmat;

%---------------Calculating fft of the stack & BPN response------
C = abs(fft(Channel'));
C(1,:) = [];
KX = C(1:block/2,:);
MI = min(KX);
MA = max(KX);
for j = 1:5,
KX(:,j) = (KX(:,j) - MI(j))/(MA(j) - MI(j));
end;
a1 = simuff(KX(1:len,:),W1,B1,trf(1,:),W2,B2,trf(2,:));
a3 = round(a1);
%---------Calculating Pearson correlation ---------------
X =Channel(Corr1(1),:)';
Y =Channel(Corr1(2),:)';
n1 = size(X);
N1 = n1(1) * sum (X.*Y);
N2 = sum(X) * sum(Y);
D1 = sqrt ( n1(1) * sum(X.^2) - sum(X)^2 );
D2 = sqrt ( n1(1) * sum(Y.^2) - sum(Y)^2 );
R1 = (N1-N2) / (D1*D2) ;

X =Channel(Corr2(1),:)';
Y =Channel(Corr2(2),:)';
n1 = size(X);
N1 = n1(1) * sum (X.*Y);
N2 = sum(X) * sum(Y);
D1 = sqrt ( n1(1) * sum(X.^2) - sum(X)^2 );
D2 = sqrt ( n1(1) * sum(Y.^2) - sum(Y)^2 );
R2 = (N1-N2) / (D1*D2) ;

V = (AMP(1)/AMP(4))/(AMP(1)/AMP(4) + AMP(2)/AMP(3));
Y = (AMP(2)/AMP(3))/(AMP(1)/AMP(4) + AMP(2)/AMP(3));

TR = [V Y R1 R2 a1(1) a1(2) a1(3) a1(4)+fact a1(5)];
Datap(i) = simuff(TR',W11,B11,trf(1,:),W21,B21,trf(2,:));

clear Channel;
clear C;
clear MI;
clear MA;
clear KX;

st = en + 1;
en = (i+1)*block;
end;


C = sort(Datap);
MinT = C(ceil(length*ProcDef.annthreshold));



j=1;
SizeRatio = ProcDef.block/block;


if SizeRatio >= 1,
   for i = 1:ProcDef.len,
      if Datap(j:j+(SizeRatio-1)) < MinT,
         ProcDef.SelStacks(i)= 0;
      else,
         ProcDef.SelStacks(i)= 1;
      end;
      j=j+SizeRatio;
   end;
end;

   
   
fprintf('\nselected %d samples\n',j);
